{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Taylor integration of the Kepler problem\n",
    "\n",
    "Here, we try to reproduce __exactly__ the [Kepler problem integration example](http://nbviewer.jupyter.org/github/JuliaDiff/TaylorSeries.jl/blob/master/examples/1-KeplerProblem.ipynb) made by Luis Benet in [JuliaDiff/TaylorSeries.jl](https://github.com/JuliaDiff/TaylorSeries.jl).\n",
    "\n",
    "The Kepler problem is the basis of planetary motion; it describes the motion of a secondary body (e.g., a planet, asteroid, comet, etc.), around a primary body (e.g., the Sun). In cartesian coordinates over the orbital plane, the Hamiltonian for the Kepler problem reads:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "H_{\\mathrm{Kepler}} &= \\frac{1}{2\\mu}(p_x^2+p_y^2)-\\frac{\\mu}{\\sqrt{x^2+y^2}}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where $\\mu=G(m_1+m_2)$, $G$ is the gravitational constant, $m_1$ is the mass of the primary body and $m_2$ is the mass of the secondary body. If we write $\\vec r_1 = (x_1,y_1)$ and $\\vec r_2 = (x_2,y_2)$, respectively, for the position of the primary and secondary body, then the vector $\\vec r = \\vec r_2-\\vec r_1$, with coordinates $\\vec r = (x, y)=(x_2-x_1,y_2-y_1)$, represents the position of the secondary body relative to the primary body; i.e., $\\vec r$ represents the so-called relative coordinates. In terms of the vector $\\vec{r}$, the position $\\vec{r}=0$ corresponds to the position of the primary body.\n",
    "\n",
    "Using Hamilton equations, we can obtain the equations of motion for the Kepler problem:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\dot x &= u \\\\\n",
    "\\dot y &= v \\\\\n",
    "\\dot u &= -\\frac{\\mu x}{(x^2+y^2)^{3/2}}\\\\\n",
    "\\dot v &= -\\frac{\\mu y}{(x^2+y^2)^{3/2}}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Note that the canonical momenta are $p_x$ and $p_y$, while $u$ and $v$ are, respectively, the  $x$ and $y$ components of the velocity; i.e., $p_x = \\mu u$ and $p_y = \\mu v$. \n",
    "\n",
    "\n",
    "First of all, we shall include all relevant packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Fontconfig warning: ignoring UTF-8: not a valid region tag\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using TaylorIntegration, Plots, LaTeXStrings\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some parameters necessary for the integration:\n",
    "\n",
    "+ $\\mu$: the gravitational parameter\n",
    "+ `q0`: the initial condition (we will select an initial condition which corresponds to elliptical motion)\n",
    "+ `order`: the order of the Taylor expansion\n",
    "+ `t_max`: the final time of the integration\n",
    "+ `abs_tol`: the absolute tolerance\n",
    "+ `n_iter`: the number of time-steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "500000"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "const μ = 1.0\n",
    "const q0 = [0.19999999999999996, 0.0, 0.0, 3.0] # a initial condition for elliptical motion\n",
    "const order = 28\n",
    "const t0 = 0.0\n",
    "const t_max = 10000*(2π) # we are just taking a wild guess about the period ;)\n",
    "const abs_tol = 1.0E-20\n",
    "const steps = 500000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we write down the equations of motion into a `function`, which here we will name `kepler!`. `q` represents the system state; i.e., the set of values of the dynamical variables at a given instant; `dq` represents the time derivatives of the components of `q`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "kepler! (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#an auxiliary array which helps with optimization:\n",
    "const r_p3d2 = Array{TaylorSeries.Taylor1{Float64}}(1)\n",
    "\n",
    "#the equations of motion for the Kepler problem:\n",
    "function kepler!(t, q, dq)\n",
    "    r_p3d2[1] = (q[1]^2+q[2]^2)^(3/2)\n",
    "    \n",
    "    dq[1] = q[3]\n",
    "    dq[2] = q[4]\n",
    "    dq[3] = -μ*q[1]/r_p3d2[1]\n",
    "    dq[4] = -μ*q[2]/r_p3d2[1]\n",
    "    \n",
    "    nothing\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Taylor integration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 43.060857 seconds (464.94 M allocations: 39.523 GB, 12.87% gc time)\n"
     ]
    }
   ],
   "source": [
    "t, q = taylorinteg(kepler!, q0, t0, 0.01, order, abs_tol, maxsteps=2); #warm-up lap\n",
    "@time t, q = taylorinteg(kepler!, q0, t0, t_max, order, abs_tol, maxsteps=steps);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10000.0,[0.2,2.72854e-9,-2.27378e-8,3.0])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t[end]/(2pi), q[end,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's extract the values of $x$, $y$, $u$ and $v$ for each time-step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, u, v = view(q,:,1), view(q,:,2), view(q,:,3), view(q,:,4);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orbital motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial conditions we selected correspond to a elliptical orbit with a relatively high eccentricity: $e=0.8$ for the initial condition `q0=[0.19999999999999996, 0.0, 0.0, 3.0]`. How does the motion of the planet/asteroid/comet looks like? Well, let's plot its orbit over the $x-y$ plane (the orbit is shown in blue; the position of the primary body is shown as a yellow dot):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scatter(\n",
    "\n",
    "[0.0], [0.0],\n",
    "label = L\"\\mathrm{Center\\; of\\; attraction}\",\n",
    "ms=10\n",
    "\n",
    ")\n",
    "scatter!(\n",
    "\n",
    "x[1:75:end], y[1:75:end],\n",
    "title = L\"\\mathrm{Orbital\\;\\;motion}\", \n",
    "xaxis = (L\"x\", (-2.0, 0.5), -2.0:0.5:2.0), \n",
    "yaxis = (L\"y\", (-0.8, 0.8), -2.0:0.5:2.0),\n",
    "label = L\"\\mathrm{Keplerian\\; ellipse}\",\n",
    "ms=0.5\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy conservation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we write the energy function of the Kepler problem; i.e., the Hamiltonian $H_{\\mathrm{Kepler}}$ in terms of $x$, $y$, $u$ and $v$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "H (generic function with 1 method)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H(x_, y_, u_, v_) = 0.5*(u_^2+v_^2)-μ/sqrt(x_^2+y_^2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, using the function `H` defined above, we calculate the energy during each time-step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = H.(x, y, u, v);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$E_0=H_{\\mathrm{Kepler}}(x_0,y_0,u_0,v_0)$, the initial value of the energy, is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.5000000000000009"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E0=E[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define $\\delta E$ as the relative error in the energy; i.e.:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\delta E(t) &= \\frac{E(t)-E_0}{E_0}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "436153-element Array{Float64,1}:\n",
       " -0.0        \n",
       "  1.77636e-15\n",
       "  1.77636e-15\n",
       "  8.88178e-16\n",
       "  8.88178e-16\n",
       "  2.66454e-15\n",
       "  3.55271e-15\n",
       "  2.66454e-15\n",
       "  3.55271e-15\n",
       "  2.66454e-15\n",
       "  3.9968e-15 \n",
       "  2.66454e-15\n",
       "  3.10862e-15\n",
       "  ⋮          \n",
       " -7.12763e-14\n",
       " -7.19425e-14\n",
       " -7.14984e-14\n",
       " -7.10543e-14\n",
       " -7.19425e-14\n",
       " -7.10543e-14\n",
       " -7.10543e-14\n",
       " -6.92779e-14\n",
       " -7.01661e-14\n",
       " -6.92779e-14\n",
       " -6.92779e-14\n",
       " -7.10543e-14"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "δE = (E-E0)/(E0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytical solution preserves the energy; i.e., $\\delta E(t)=0$ for the analytical solution. Thus, we expect our solution to be close to zero. So, even if the $\\delta E$ is not perfectly zero during the whole integration, it is comparable to Julia's `Float64` machine-epsilon.\n",
    "\n",
    "Below, we plot $\\delta E$ in units of Julia's `Float64` machine-epsilon as a function of time, $t$. In Julia, the machine-epsilon for the `Float64` type has a value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.220446049250313e-16"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eps(Float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(\n",
    "\n",
    "t/(2π), δE/eps(Float64),\n",
    "title = L\"\\mathrm{Energy\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\", \n",
    "xaxis = (L\"t\\mathrm{\\,(orbital\\;\\;periods)}\"),\n",
    "yaxis = (L\"\\delta E \\mathrm{\\;\\;(machine\\;\\;epsilons)}\"),\n",
    "label = L\"\\mathrm{Energy\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\"\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, how does the energy error distribute around zero?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[34mINFO: binning = 20\n",
      "\u001b[0m"
     ]
    },
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "histogram(\n",
    "\n",
    "δE/eps(Float64),\n",
    "title = L\"\\mathrm{Distribution\\;\\;of\\;\\;energy\\;\\;relative\\;\\;error}\", \n",
    "xaxis = (L\"\\delta E\"),\n",
    "yaxis = (L\"\\mathrm{Number\\;\\;of\\;\\;}\\delta E\\mathrm{\\;\\;values\\;\\;within\\;\\;bin\\;\\;range}\"),\n",
    "leg = false,\n",
    "nbins=20\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy error behaves roughly as a random walk, which means that the numerical error in the energy is dominated by rounding errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Angular momentum conservation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analogously to the energy, we now focus on the angular momentum, which is preserved by the analytical solution too. The value of the angular momentum is given by\n",
    "$$\n",
    "\\begin{align}\n",
    "L &= xv-yu\n",
    "\\end{align}\n",
    "$$\n",
    "We write the angular momentum function as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ang_mom (generic function with 1 method)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ang_mom(x_, y_, u_, v_) = x_.*v_-y.*u_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the angular momentum during each time-step is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = ang_mom(x, y, u, v);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$L_0=x_0v_0-y_0u_0$, the initial value of the angular momentum, is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5999999999999999"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L0 = L[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define $\\delta L$ as the relative error in the angular momentum; i.e.:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\delta L(t) &= \\frac{L(t)-L_0}{L_0}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "δL = (L-L0)/L0;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just as we did for the energy, we now plot $\\delta L$ in units of `eps(Float64)` vs $t$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(\n",
    "\n",
    "t/(2π), δL/eps(Float64),\n",
    "title = L\"\\mathrm{Angular\\;\\;momentum\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\", \n",
    "xaxis = (L\"t\\mathrm{\\,(orbital\\;\\;periods)}\"),\n",
    "yaxis = (L\"\\delta L \\mathrm{\\;\\;(machine\\;\\;epsilons)}\"),\n",
    "label = L\"\\mathrm{Angular\\;\\;momentum\\;\\; relative\\;\\; error\\;\\; vs\\;\\; time}\",\n",
    "leg=false,\n",
    "color=:orange\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we see that the angular momentum relative error is comparable to `eps(Float64)`; the maximum variation is ~$150$ machine epsilons. How does the distribution of this error look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[34mINFO: binning = 20\n",
      "\u001b[0m"
     ]
    },
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "histogram(\n",
    "\n",
    "δL/eps(Float64),\n",
    "title = L\"\\mathrm{Distribution\\;\\;of\\;\\;angular\\;\\;momentum\\;\\;relative\\;\\;error}\", \n",
    "xaxis = (L\"\\delta L\"),\n",
    "yaxis = (L\"\\mathrm{Number\\;\\;of\\;\\;}\\delta L\\mathrm{\\;\\;values\\;\\;within\\;\\;bin\\;\\;range}\"),\n",
    "leg = false,\n",
    "nbins=20,\n",
    "color=:orange\n",
    "\n",
    ")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The distribution of $\\delta L$ shows a peak near zero and roughly symmetrical around this value. This means that the error in the angular momentum is dominated also by rounding errors of floating-point arithmetic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we reproduce __exactly__ the last plot shown in the original version of this example, authored by Luis Benet and included in `JuliaDiff/TaylorSeries.jl`'s [Kepler problem integration example](http://nbviewer.jupyter.org/github/JuliaDiff/TaylorSeries.jl/blob/master/examples/1-KeplerProblem.ipynb) jupyter notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A $\\delta E$, $\\delta L$ plot vs $t$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(\n",
    "\n",
    "[t/(2π), t/(2π)], [δE/eps(Float64), δL/eps(Float64)],\n",
    "title = L\"\\delta E\\mathrm{\\;\\;(blue),\\;\\;}\\delta L\\mathrm{\\;\\;(green)\\;\\;vs\\;\\;time}\", \n",
    "xaxis = (L\"t\\mathrm{\\,(orbital\\;\\;periods)}\"),\n",
    "yaxis = (L\"\\delta E\\mathrm{,\\;\\;}\\delta L\\;\\;\\mathrm{(machine\\;\\;epsilons)}\"),\n",
    "leg=false,\n",
    "\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.3-pre",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
